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Abstract 


) 

A  nonlinear,  Square-Root  Variable-Metric  optimization 
technique  is  used  to  extract  time  histories  for  the  lift 
coefficient,  drag  coefficient,  and  bank  angle  control 
variables  from  radar  data  of  the  Space  Shuttle  Test  Flight. 

I.  This  optimization  technique  solves  for  these  control 
variable  histories  by  minimizing  a  weighted  least  squares 
performance  index  of  the  calculated  versus  experimental 
state  equation  time  histories.  The  state  model  is  a  three 
degree  of  freedom  representation  of  aircraft  motion. 

Three  intervals  of  data  are  evaluated.  In  each  inter¬ 
val,  polynomial  functions  represent  the  control  variables, 
and  the  peformance  index  is  dependent  only  on  the  values 
of  the  polynomial  coefficients.  The  larger  of  the  first 
two  intervals  contained  the  shorter  interval  so  that  the 
effect  of  changing  polynomial  order  and  interval  length 
could  be  explored.  The  third  interval  compared  the  calcu¬ 
lated  parameter  histories  with  NASA's  BET  parameter  histories 
with  good  correspondence.  For  all  three  intervals  the 
calculated  state  histories  closely  matched  the  experimental 
state  histories. 

These  results  show  that  given  an  appropriate  state 
model,  this  optimization  technique  is  useful  for  quickly 
obtaining  good  and  inexpensive  estimates  of  certain  desired 
parameters. 

Approved  for  public  release ;  distribution  unlimite 
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NONLINEAR  PARAMETER  IDENTIFICATION: 

LIFT  COEFFICIENT,  DRAG  COEFFICIENT,  AND 
BANK  ANGLE  HISTORIES  FOR  THE  SPACE  SHUTTLE 
TEST  FLIGHT  I 

I  Introduction 

Evaluating  the  performance  of  a  new  and  technologic¬ 
ally  advanced  engineering  system  requires  a  thorough,  and 
often  redundant,  study  of  test  data.  Agreement  between 
these  different  studies  helps  to  reinforce  any  conclusions 
derived  from  the  tests.  The  Space  Shuttle  Test  Flights 
will,  accordingly,  be  very  well  studied  and  documented. 

The  Space  Shuttle  trajectory  is  computed  several  different 
ways,  and  this  data  enables  important  parameters,  such  as 
lift  coefficient,  drag  coefficient,  and  bank  angle  histories 
to  be  determined.  These  parameter  histories  give  insight 
into  the  performance  of  the  Shuttle,  so  that  expected 
performance  parameters  can  be  evaluated  for  their  accuracy. 

NASA  computes  a  "Best  Estimated  Trajectory"  (BET) 

(Ref  7),  containing  estimates  of  the  parameter  histories, 
by  using  the  data  from  the  Aerodynamic  Coefficient 
Instrumentation  Package  (ACIP) ,  which  is  on  board  the  Space 
Shuttle,  and  various  radar  data.  The  ACIP  gives  gyro  and 
accelerometer  data,  whereas  the  radar  gives  optical  and 
telemetry  data.  This  data  is  combined,  through  a  linear 


stochastic  estimation  process,  to  produce  an  optimal  six 
degree  of  freedom  (DOF)  model  of  the  Space  Shuttle's  trajec¬ 
tory  (Ref  3,  4,  5,  6). 

Although  NASA's  BET  is  a  good  method  for  estimating 
the  trajectory  and  motion  of  the  Space  Shuttle,  it  may  take 
a  few  months  to  synthesize  all  of  the  available  data  into  a 
BET.  Yet  the  radar  data,  such  as  the  data  from  the  Multi¬ 
sensor  Space  Position  Report  (Ref  8:1),  is  sufficient  to 
calculate  three  DOF  trajectories  and  performance  parameters, 
such  as  Cl,  Cd,  and  /*-  .  This  report  proposes  a  simple  and 
inexpensive  computer  algorithm  to  evaluate  this  data  and 
come  up  with  a  good  estimate  of  some  of  the  performance 
parameters  of  the  Space  Shuttle.  This  could  provide  an 
on-line  estimate,  which  evaluates  the  radar  data  as  it  comes 
in.  Its  results  could  be  compared  to  the  BET  in  post  data 
analysis,  which  would  help  locate  possible  software  errors 
in  the  two  estimation  procedures  and  be  used  to  confirm 
wind  tunnel  data.  The  immediate  results  would  also  be 
useful  in  evaluating  the  trajectory  should  anything  go 
wrong  during  the  final  landing  phase  of  reentry. 

Problem  Statement  and  Method 

This  report  tests  a  parameter  estimation  method  which 
can  be  used  to  evaluate  the  Space  Shuttle  radar  data  as  it 
becomes  available.  The  method  proposed  is  a  slightly 
modified  Davidon-Fletcher-Powell  (DFP)  optimization  technique 
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derived  by  Williamson  and  Hull  (Ref  4:1).  The  dynamical 
system  from  which  the  parameters  are  identified  can  be 
closely  approximated  without  the  need  for  linear  stochastic 
noise  modeling  by  using  a  deterministic  nonlinear  state 
model  (Ref  4) .  The  DFP  variable-metric  method,  using  such 
a  model,  is  an  efficient,  quasi  second-order  numerical 
algorithm  that  minimizes  a  performance  index.  This  perfor¬ 
mance  index  is  calculated  from  a  weighted  least  squares  fit 
of  the  calculated  onto  the  experimental  (radar  data)  state 
histories.  The  calculated  state  histories  are  dependant 
upon  the  parameter  histories  through  the  coefficients  of  a 
polynomial  fitted  to  each  parameter  history.  These  coeffi¬ 
cients  are  the  variables  that  are  changed  to  minimize  the 
performance  index.  The  parameter  histories  that  this  report 
identifies  are  the  lift  coefficient,  the  drag  coefficient, 
and  the  bank  angle  histories. 

The  data  used  to  calculate  the  performance  index 
came  from  two  sources.  The  first  set  of  data  came  from  the 
Mulitsensor  Space  Position  Report  (Ref  8  -  1 ) .  Its  trajec¬ 
tory  histories  supply  sufficient  information  to  derive  the 
time  histories  of  the  parameters  CL,  CD,  and  through  some 
optimization  technique,  such  as  the  algorithm  in  this  report. 
As  this  algorithm  was  being  evaluated  using  the  Multisensor 
Report  data,  the  NASA  BET  data  became  available.  Although 
the  BET  data  was  itself  derived  from  a  parameter  estimation 
process,  it  contains  trajectory  history  data  as  in  the 
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Multisensor  Report.  The  parameter  histories  obtained  by 
the  algorithm  using  these  trajectory  histories  were  compared 
to  the  BET  parameter  histories.  The  degree  of  correspondence 
between  these  parameter  histories  provides  a  measure  of  the 
accuracy  of  the  algorithm. 

Organization 

Chapter  II  discusses  the  procedure  used  to  extract  the 
parameters  from  the  radar  data.  This  will  include  a  descrip¬ 
tion  of  the  Square -Root  Variable-Metric  algorithm  and  an 
explanation  of  the  mating  of  the  state  equations  and  data. 
Chapter  III  explains  the  problems  and  pitfalls  of  manipulating 
the  algorithm  to  provide  the  best  results.  Chapter  IV  presents 
the  results  and  interprets  their  significance.  Chapter  V 
summarizes  the  -report  and  gives  an  explanation  of  the  useful¬ 
ness  of  this  parameter  identification  method. 


II  Procedure 


Methods  of  Parameter  Identification 

Many  techniques  have  evolved  that  solve  for  variables 
of  interest  from  some  collection  of  input.  Each  method  has 
proven  useful  at  some  time,  although  the  degree  of  accuracy 
between  methods  varies  greatly.  For  instance,  the  BET  gives 
accurate  parameter  estimations  given  all  of  the  available 
data.  Yet  it  was  a  time-consuming  process  to  get  the  BET. 

On  the  other  hand,  there  are  gradient  optimization  techniques 
that  are  quick  and  inexpensive,  yet  only  crudely  approximate 
the  parameters.  But  there  are  certain  common  characteristics 
of  parameter  estimation  methods  for  a  dynamical  system.  Hull 
and  Williamson  (Ref  4:1)  explain: 

In  general,  the  parameter  estimation  problem 
for  a  dynamical  system  consists  of  finding  the 
values  of  a  set  of  p  constants  C  and  n  initial 
states  XQ  (assuming  to=0)  such  that  the  solution 
of  the  dynamical  system 


F  (i  ,X  C) 


fits  the  experimental  data.  The  data  is  fit  by 


minimizing  the  weighted  least  squares  performance 
index 


0  =  I  (  X  -  &S  W  (  Y*  -  Cr  J 

K 


(2) 


were  the  subscript  k  denotes  the  time  at  which 
the  experimental  data  is  taken  and  W  is  a  positive- 
definite  diagonal  matrix  whose  elements  are  the 
inverses  of  the  measurement  convariance.  The 
quantity  Y  is  a  q-vector  of  the  experimental  values, 
and  G  is  q-vector  of  computed  values  which  is  related 
to  the  state  of  the  dynamical  system,  that  is, 

G  =  G(t,X) . 


The  technique  of  parameter  identification  proposed  in 
this  report  is  a  compromise  between  the  crude  and  the  sophis¬ 
ticated  methods.  It  is  a  nonlinear  method,  which  implies 
that  the  state  equations  would  exactly  map  the  motion  if  the 
three  DOF  model  was  an  exact  representation  of  the  real 
world.  The  Square  Root  Variable-Metric  algorithm  of 
Williamson's  (Ref  9:107)  chooses  [x0  C]T,  where  XQ  is 

the  initial  state  vector  and  C  is  a  vector  of  constants 
associated  with  the  parameter  histories,  to  minimize  the 


weighted  least  squares  performance  index-.  The  algorithm  is: 


1.  Guess  £  and  a  positive-definite,  symmetric 
matrix  S.  Compute  J((3),  J|»(p)  (the  derivative 
of  the  performance  index  with  respect  to  a ) 
and  H  =  SST. 

2.  Compute  the  change  in  ^  from  the  formula 

A p>  =  -  -c  HJT  (3) 

using  a  one-dimensional  search  to  determine 
the  scalar  step  size  ®<  which  minimizes  J  in 
the  -HJt  direction.  This  search  yields  a  new 
function  value  J,  where 

J  =  J(  (i  +  a/S)  (4) 

3.  Compute  a  new  gradient  and  form  the  difference 

^jT  =  j?  ‘  (5) 

4.  Compute  a  new  metric  H  from  the  relation 

H  =  SST  (6) 

where 

S  =  S [I  -  (ty/A)STYYTS]  (7) 

A  =  YT  H.f  Y  =KjT  + 

■V  =  [l+d-B/A)'^]  /  (B/A)  ,  B  =  YTSSTY 

If  1)  is  imaginary,  use  ~u  =  1  in  the  above 
expression  for  S. 

5.  If  the  process  has  not  converged,  repeat  the 
procedure.  Convergence  has  occurred  when  cx— *1, 

A — *0  and  B  — *  0  simultaneously. 


Since  the  metric  H  is  calculated  as  the  S  matrix  times 
its  transpose,  it  will  always  be  positive  definite.  The 
positive  definiteness  of  H  means  the  performance  index  will 
decrease  on  every  iteration,  regardless  of  the  accuracy  in 
calculating  the  S  matrix.  Likewise,  the  one-dimensional 
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alpha  search  need  not  be  extremely  accurate,  as  is  necessary 
with  the  unmodified  DFP  algorithm,  to  maintain  a  positive 
definite  H. 


Numerical  Integration  and  Differentiation 

The  method  used  to  numerically  calculate  J(|g)  is  as 
follows : 

1.  Guess  values  for  XQ  and  C  to  give  p  =  [J^C]T. 

2.  Compute  Y0  -  G(t0,  XQ)  from  guessed  Xq.  (Note 
that  in  this  report  the  initial  state  vector 
XQ  was  assumed  to  be  equal  to  Y0  and  &  was 
therefore  -  [C]T.) 

• 

3.  Integrate  X  =  F(t,  X,  C)  from  tQ  to  t^ .  This 
enables  Y.  -  G(ti,,XiL)  to  be  computed.  In  this 
manner,  integrate  X  and  calculate  Y^  -  G(tv,  XK) 
until  the  final  time  (the  last  observation)  is 
reached. 

4.  Using  Eq.  (2),  compute  J((3). 


The  gradient  of  J(p),  (p ) ,  is  a  vector  calculated 
through  numerical  differentiation  by  perturbing  each  element 
of  £  and  differencing  the  performance  indexes  that  are  calcu¬ 
lated  using  the  perturbed  (h  vectors.  The  central  difference 
formula 


\ --  [J (f  I,  fix,  •  ■ • ' *  >  •  •  *  >  Pr.n)  -  JO.,  ^  , . . . ,  j]  /( 2  Wj  ) 


(8) 


is  the  differencing  technique  used  in  this  report,  where 
is  the  perturbation  in  the  element  of  .  This  means 

that  each  element  of  requires  two  calculations,  and  hence 
two  integration  loops,  for  J.  The  central  difference  formula 
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mmm 


w 


r 

[■  ,  was  chosen  because  it  is  uncomplicated  yet  provides  a 

|  truncation  error  of  order  ^ (*>*"  .  Therefore,  the  smaller  the 

i 

j>  the  better,  but  there  is  also  a  limit  to  how  small 

can  be  due  to  computer  round-off  error. 


[  Hull  and  Williamson  (Ref  4)  say  that  "in  order  to  have 

> 

'  the  greatest  number  of  correct  significant  figures  in  a 

derivative,  it  is  necessary  to  select  the  C3  which  causes 

5 

the  largest  change  in  J  while  keeping  the  truncation  error 
of  the  derivative  in  check."  J(  p-  ±  $  (b-  ),  and  hence  , 
l  are  computed  using 


5  ft  -  ^  l  ft  I  if  I  ft,  I  >  1  (9) 

6  ft  =  €  if  l  ft  U  1 


Fig.  1.  Flat-Earth  Coordinate  System 

These  equations  are  derived  from  a  flat-Earth  coordinate 
system  (Fig.  1)  where  x3  is  the  altitude,  x-^  and  x2  are  the 
planar  coordinates  orthogonal  to  the  altitude  such  that  the 
coordinate  system  is  right  handed.  The  directions  of  x-^  and 
x2  depend  on  which  data  set  is  used.  The  state  X4  is  the 
magnitude  of  the  velocity,  X5  is  the  azimuth,  and  xg  is  the 
flight  path  angle. 

Data  Manipulation 

The  algorithm  requires  an  experimental  trajectory  (YK) 
as  an  input.  The  YK  j.s  a  time  history  of  the  states.  The 
radar  data  can  be  manipulated  to  give  state  vector  values 
at  each  time  increment  k.  This  report  used  two  sources  of 
radar  data  to  yield  YK,  the  Multisensor  Report  data  and  the 
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BET  data .  The  data  from  the  Mulitsensor  Report  (Ref  8) 
was  with  respect  to  an  Earth  Surface  Fixed  (ESF)  Cartesian 


~  North 


A  1 1  rv-s  t  V\ 


X,D-  Fait 


Fig.  2.  Earth  Surface  Fixed  Coordinate  System 


coordinate  system  (Fig.  2)  which  rotates  with  the  Earth. 

Because  the  time  intervals  evaluated  were  very  short  (less 

than  two  minutes  in  every  case) ,  the  rotation  of  the  Earth 

was  ignored  in  transforming  the  radar  data  into  Y  .  As 

K 

shown  in  Figure  2,  the  ESF  coordinate  system  had  the 
beginning  of  runway  23,  Edwards AFB,  as  the  origin,  with  the 
xid  axis  positive  East,  the  x2D  axis  positive  North,  and 
the  x3D  axis  positive  up.  For  YK ,  the  states  x1  and  x2 
were  given  by  xlD  and  x2D  data,  respectively.  The  azimuth, 
X5D'  was  measured  from  North  with  the  clockwise  direction 


being  positive,  and  the  flight  path  angle,  xgD,  was  the 
angle  between  the  velocity  vector  and  the  (x1D,  x2D)  plane. 
For  Yr  the  states  x5  and  xg  were  given  by  x5D  and  xgD, 
respectively,  and  x4  was  given  by  the  magnitude  of  the 
velocity,  which  the  Multisensor  Report  also  provided. 

The  only  state  that  needed  to  be  calculated  from  the 


Fig.  3.  Runway  Coordinate  System 


given  data  was  x^.  The  Multi sensor  Report  also  provided  the 
altitude  of  the  Space  Shuttle.  Therefore  X3  was  given  by 
the  Space  Shuttle  altitude  minus  the  altitude  of  the  runway. 
The  BET  data  uses  a  similar  coordinate  system  (Fig.  3) . 
The  x^d  axis  is  positive  in  the  direction  of  the  runway;  the 
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positive  X2D  axis  makes  up  the  third  coordinate  for  the  right 
handed  coordinate  system  with  the  x^  axis  positive  down  (Ref  1)  . 
Again,  x^,  (asimuth)  is  measured  from  North,  clockwise 
positive,  and  XgD  was  the  angle  between  the  velocity  vector 
and  the  (xiD,  X2D)  plane.  To  arrive  at  state  elements  xj_  and 
X2  of  Yk  a  rotation  of  25.5  degrees  was  added  to  the  x^  term 
in  each  of  the  x^  and  &2  equations.  This  gave  the  xj  and 
X2  equations  a  corrected  azimuth  measured  from  the  X2D  axis. 

Although  the  Multisensor  Report  provided  the  altitude 
of  the  Space  Shuttle,  for  the  BET  data  the  altitude,  and 
therefore  X3,  the  altitude  of  the  Space  Shuttle  measured 
from  the  origin  at  runway  23,  needed  to  be  calculated  using 
the  cosine  formula  for  triangles 


Figure  4  shows  the  geometry  that  this  equation  was  derived 
from.  This  equation  is  applicable  as  long  as  X3D  is  negative 
(Space  Shuttle  is  above  the  horizon) .  Since  this  report  is 
concerned  with  data  from  radar  stations  along  the  VTest  Coast, 
the  problem  of  X3J3  positive  was  not  considered.  The  X5  and 
x6  elements  of  YK  are  given  by  x5D  and  x6D  respectively. 
Finally,  the  magnitude  of  the  velocity,  X4 ,  is  calculated 
as  the  square  root  of  the  sum  of  the  squares  for  x^p*  *2D' 
and  X3D  (which  are  provided  by  the  BET) . 
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Parameters 


The  parameter  histories  of  interest  are  C L,  CD,  and  . 
The  algorithm  requires  initial  guesses  for  these  parameter 
histories.  To  solve  for  the  initial  time  history  for  CT  the 
equation  of  motion  in  the  vertical  direction 

l_  c 05  (V)  COS  (V)  laJ~  ~  fh  Vt  (17) 

was  used.  The  cos^<)  term  was  chosen  to  be  equal  to  one. 

This  approximation  is  based  upon  the  knowledge  that  ^  for  a 
heavy  glider  like  the  Space  Shuttle  should  not  exceed  45 
degrees,  and  in  most  cases  it  will  be  below  30  degrees. 
Therefore,  at  worst,  this  is  a  70%  accurate  approximation. 
Furthermore,  vz(t-),  which  is  the  acceleration  in  the  Z  (up) 
direction  at  time  t^ ,  can  be  approximated  by  the  central 
difference  formula 

vt(o=  (v,(u  -  - 1.. )  ns) 

Therefore,  representing  L  as  <ip^«)v\,  5  CjC< .  '  we  have 

C(Q  «  ~  vjt.,  )/($(<>*,  ~  t-.))  +  2-~]  t  (19) 

fU. ,  s  j  c o%( ,) 

where  n  is  the  length  of  the  time  interval  of  interest 
divided  by  the  step  size  taken  between  data  points  (n  is  the 
total  number  of  data  points  for  each  of  CL,  Cp,  and  jla)  and 
where  Y(tj)  is  the  flight  path  angle  at  time  t( .  The  air 
density  (p)  history  is  also  derived  from  a  Least  Squares 
polynomial  fit  to  the  air  density  data  given  by  References 
7  and  8 . 


... 
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For  the  Multisensor  Report  data,  the  CD  time  history 
guess  was  derived  from  lift  over  drag  graphs  in  Reference  1. 
The  CD  history  for  the  BET  data  was  calculated  using  the  lift 
over  drag  (L/D)  data  already  available  from  the  BET.  However 
the  initial  Cl  and  Cd  histories  are  calculated,  the  closer 
they  are  to  the  true  CL  and  CD  the  faster  the  algorithm  will 
converge.  These  initial  histories  will  be  compared  to  the 
converged  histories  to  determine  whether  or  not  they  provide 
a  good  initial  guess  for  the  estimation  method. 

For  both  sets  of  data  the  ju.  initial  time  history  was 
derived  from  the  Xg  state  equation,  again  using  a  discrete 
time  central  difference  to  represent  the  derivative  such 
that 


/*(t,  )= 


UJ 


tc4x  v  < 


>(X«> 


'  s  v  C ) 


(20) 


Note  that  Y,X  ,  and  v  represent  flight  path  angle,  azimuth, 
and  magnitude  of  the  velocity,  respectively,  as  did  x6 f  x5 , 
and  X4  in  the  X5  equation. 

Presently,  ^  is  made  up  of  three  [nxl]  vectors  CL,  CD, 
and  yLK.  (as  mentioned  previously,  n  is  the  number  of  data 
points  for  each  of  CL,  CD,  and  v),  hence,  is  a  r3nxl] 
vector.  Since  the  performance  index  is  numerically  differ¬ 
entiated  with  respect  to  each  element  of  ji  ,  a  smaller 
number  of  elements  of  ^  will  greatly  reduce  computation 
time.  Because  n  is  the  time  interval  length  divided  by  the 
time  increment  size,  if  n  is  small  then  either  the  increment 
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size  is  large  or  the  interval  length  is  short.  Since  it 
is  desired  to  fit  the  data  points  to  a  smooth  "true"  function 
of  the  parameters,  a  relatively  short  increment  size  is  re¬ 
quired  whether  the  interval  is  long  or  short.  Therefore,  the 
size  of  n  is  predominately  dependent  upon  the  length  of  the 
interval.  To  reduce  the  computational  load  by  decreasing 
the  size  of  n  would  then  mean  a  short  interval.  This  would 
lead  to  patching  short  intervals  together  in  order  to  obtain 
good  parameter  histories  for  a  reasonable  length  of  time. 

To  make  n  as  large  as  possible  and  still  cut  down  on 
the  length  of  ,  a  Least  Squares  polynomial  fit  is  made  to 
the  Cl,  Cd,  and  y*.  histories.  Since  the  true  CL,  CD,  and  ^ 
time  histories  will  be  continuous  functions  of  time,  a 
Least  Squares  polynomial  fit  to  the  discrete  time  parameter 
data  points  is  a  reasonable  model  for  their  true  histories. 
The  resulting  coefficients  are  used  as  p  and  are  called  upon 
to  calculate  CL»  Cp,  and  ,  at  any  given  time,  for  use  in 
the  state  equations . 

A  programming  decision  must  be  made  when  deciding  the 
order  for  the  polynomials.  Reducing  the  size  of  &  without 
sacrificing  interval  length  was  the  initiative  for  using  a 
Least  Squares  polynomial  model  of  the  parameters,  therefore 
it  is  desired  to  fit  the  interval  with  the  smallest  order 
polynomial,  thereby  decreasing  the  length  of  ^  ,  and  yet 
maintain  an  accurate  functional  model  of  the  parameter 
histories.  A  rough  plot  of  the  initial  parameter  histories 
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will  give  an  idea  what  order  polynomial  will  best  approximate 
the  true  parameters.  Thus,  within  this  chapter  the  methodology 
for  identifying  parameters  is  explained  as  a  weighted  least 
squares  minimization  process.  The  Square-Root  Variable-Metric 
algorithm,  used  in  this  report,  minimizes  the  performance 
index,  J,  by  a  modified  DFP  optimization  technique  which 
is  quasi  second-order  with  respect  to  the  parameters  of 
interest  (^>)  .  The  derivative  of  J  with  respect  to  ^  is 
calculated  numerically  using  the  central  difference  formula. 

The  state  equation  model  of  the  three  DOF  motion  of  the  Space 
Shuttle  was  a  flat  Earth  approximation.  The  experimental 
state  histories  were  extracted  from  the  radar  data  according 
to  the  coordinate  system  used.  The  initial  parameter  histories 
were  derived  from  the  equation  of  motion  in  the  vertical 
direction  (for  C^) ,  Reference  1  (for  CD) ,  and  the  x5  state 
equation  (for  /a.  )  .  Finally,  the  method  of  least  squares  was 
used  to  fit  polynomials  to  the  parameter  histories,  and  the 
length  of  the  vector  was  reduced  by  using  the  coefficients 


of  these  polynomials  as 


Ill  Programming  Considerations 


After  fitting  the  appropriate  polynomials  to  the 
guessed  parameter  histories,  as  was  discussed  in  Chapter  II, 
the  algorithm  is  ready  to  find  the  optimal  .  This 

chapter  will  expand  upon  some  of  the  more  delicate  program¬ 
ming  characteristics  of  the  algorithm. 

One  Dimensional  Search 

The  one  dimensional  search  finds  the  value  of  o<  that 
minimizes  J  in  the  -HJT  direction. 


min  J  (  /g  +  A  /3>) 


where 


(21) 


A  /3  =  —  °<  HJt 

Note  that  it  is  this  -HJT  that  supplies  the  algorithm  with 
quasi  second-order  convergence  characteristics,  since  as  H 
is  updated  it  approaches  the  inverse  of  the  Hessian  matrix. 
There  are  many  methods  for  arriving  at  the  correct  cx  but 
all  methods  have  certain  inherent  problems.  The  algorithm 
becomes  very  sensitive  to  very  small  changes  in  (. 3  as 

convergence  is  approached.  Therefore,  it  is  important  to 
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keep  from  choosing  a  /3  that  give  a  value  for  J  too  far 

beyond  the  minimum  in  the  -HJ  direction  since  even  a 

small  a  p  added  to  the  can  cause  great  changes  in 

J.  This  could  result  in  improper  updating  of  the  S  matrix 
and  likewise  lead  to  searches  in  the  wrong  direction  on 
following  iterations. 

The  method  used  to  calculate  in  this  report 

(Ref  2)  involved  setting  =><  •  equal  to  c*.  +  £  ,  where 

c  *  \ 

i=l,2,...;  and  °<r0  ,  the  initial  *><  ,  is  0.,  and  S 

is  a  small  constant  of  magnitude  less  than  0.1.  This 
new  <x.  is  used  to  calculate  (2 ■  =  a  .+  a  , 

and,  hence,  J(-  (  (i-  )  •  This  J ■  (  )  is  compared  to 

Ji-1  (  /Sy. ,)  •  If  is  less  than  Jt*_i,  several 

decisions  must  be  made  before  the  minimizing  can 

be  found.  If  the  difference  between  and  J('_i  is  less 
than  some  tolerance,  then  J  can  be  considered  the 
minimum.  If  the  difference  is  greater  than  this 
tolerance,  the  size  of  £  is  decreased  and  the 
search  is  restarted  from  J*_i.  This  procedure  of  decreasing 
and  restarting  the  search  continues  until  (J^-i 
J(*)/J'-l  is  less  than  the  required  tolerance.  Note  that 
txr  must  indicate  the  sum  of  all  initialized  and  reinitial¬ 
ized  step  sizes.  Also,  since  it  will  remain  positive 
definite,  the  tolerances  need  not  be  extremely  small  to 
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maintain  good  convergence  of  the  algorithm. 


Convergence  Criteria 

The  convergence  criterion  is  based  on  (JD  -  _i)/J0 

less  than  some  tolerance,  where  Jn  and  J.  ,  are  minimum 

u  t  -1 

values  of  the  performance  index  from  two  successive  search 
directions.  If  (JQ  -  Jj-i)/J0  less  than  the  prescribed 
tolerance,  and,  likewise,  A  — *  0,  B — *  0,  and  <x  — >1,  then 
the  algorithm  has  converged.  As  convergence  is  approached, 
the  minimum  J,  J<‘_i,  may  be  the  same  as  the  value  of  J0  if 
the  one  dimensional  search  only  takes  one  step  before  its 
tolerance  is  satisfied.  For  this  case,  if  the  percentage 
difference  (JQ  -  is  greater  than  the  prescribed 

tolerance,  then  this  tolerance,  and  the  tolerance  for  the 
one  dimensional  search  are  decreased  and  the  one-dimensional 
search  is  restarted.  If  the  percentage  difference  is  less 
than  the  prescribed  tolerance,  and  A,  B,  and  have  not 

converged,  then  the  S  matrix  is  reinitialized  to  the  identity 
matrix,  I,  and  a  gradient  step  is  taken.  This  is  done  on 
the  chance  that  H  is  not  converging  to  the  inverse  Hessian 
matrix  and  needs  to  be  reinitialized.  Generally,  the 
algorithm  has  reached  convergence  before  this  step  is  used. 

Differentiation 

Another  programming  consideration  is  the  choice  of  the 
best  6  for  the  numerical  di fferention  of  J ( R  ) .  Hull 


22 


and  Williamson  (Ref  4)  explain  a  method  for  calculating  the 
best  6  which  involves  taking  second  numerical  derivatives. 
Instead,  to  simplify  the  algorithm  further,  derivatives 
were  calculated  with  several  different  values  of  .  A 

value  of  C  r  \,0£~H  gave  the  largest  number  of  significant 
digits  in  the  derivative  without  it  being  affected  by 
computer  roundoff,  making  this  the  £  of  choice  for  this 
algorithm. 

Nondimens ional  Coefficients 

For  higher  order  polynomials,  such  as  the  sixth  order 
polynomial  used  to  fit  the  parameter  histories,  the  relative 
magnitude  between  coefficients  may  be  very  large.  Since 
these  coefficients  are  all  differentiated  using  the  same 
size  ,  the  elements  of  the  J  vector  may  likewise  be  of 

widely  varying  magnitudes.  This  makes  the  algorithm  more 
sensitive  to  some  elements  than  others,  thereby  retarding 
the  overall  update  of  S.  Also,  this  disparity  in  J  element 
magnitude  will  increase  the  truncation  error,  hence,  weaken¬ 
ing  the  ability  of  the  algorithm  to  converge. 

This  algorithm  therefore  used  a  nondimensional 
vector  to  make  the  magnitudes  of  the  polynomial  coefficients 
more  "even",  thereby  making  the  Jj  elements  about  the  same 
order  of  magnitude  also.  This  involved  multiplying 
the  coefficients  by  the  final  time  of  the  interval  to  the 
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power  that  the  coefficient  has  in  the  polynomial  (for  instance, 

a  second  order  polynomial,  CQ  +  C^t  +  C2t^,  would  have  CQ, 

2 

C^tf  and  C2tf  as  nondimensional  coefficients)  .  These 
elements  of  p  had  to  be  redimensionalized,  by  dividing  by 
the  appropriate  power  of  tf,  before  calculating  the  CL(t{), 
Cp(t;),  and  ^c^(tt)  values  called  for  in  the  state  equations 
during  the  integration. 

Weighting  Matrix 

The  weighting  matrix  had  to  be  adjusted  to  find  the 
best  choice  of  diagonal  elements  for  the  fastest  convergence. 
The  Multisensor  Report  gives  an  idea  of  the  uncertainties 
to  expect  of  radar  data  (Ref  8)  as  shown  in  the  Appendix. 

The  use  of  a  diagonal  weighting  matrix  is  a  restriction  to 
the  full  inverse  measurement  covariance  matrix  weighting. 
Although  the  statistical  description  of  the  radar  data  state 
errors  are  correlated,  a  diagonal,  uncorrelated,  error  model 
was  used  to  simplify  the  computations,  and  hence,  speed 
the  convergence  of  the  algorithm.  The  fastest  convergence 
was  achieved  by  weighting  the  values  of  x^  and  Xg  (angles) 
up  to  eight  orders  of  magnitude  more  than  the  x^,  x.2>  and 
x3  (distances)  states,  and  up  to  four  orders  of  magnitude 
more  than  x4  (velocity) .  Comparisons  are  made  in  Chapter 
IV  of  runs  using  different  weights. 
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IV  Results 


Two  runs  were  made,  using  Multisensor  Report  data,  to 
test  the  convergence  characteristics  of  the  algorithm.  Once 
the  algorithm  was  working  properly,  a  run  using  BET  data 
was  fully  converged  and  analyzed.  The  two  test  runs,  (call 
them  Run  1  and  Run  2),  were  straight  landing  approaches. 
Therefore  the  equation  (rate  of  change  of  azimuth)  was 
nearly  zero,  implying  very  small  bank  angles.  The  BET  data 
was  taken  from  a  turn.  This  arrangement  of  data  runs  allowed 
evaluation  of  the  algorithm  characteristics  during  a  rela¬ 
tively  simple  flight  history  before  testing  the  performance 
of  the  algorithm  during  a  period  of  active  maneuvering. 

Run  1 

The  first  run  spanned  a  seventy-eight  second  time  inter¬ 
val,  ranging  in  altitude  from  8667.5  feet  to  14.9  feet. 

Table  1  lists  the  calculated  versus  experimental  state 
histories,  and  the  difference  between  them,  element  by  element. 
This  run  used  500  seconds  of  Central  Processor  (CP)  time 
to  converge  to  A  =  9.2  E  -  6 ,  B  =  3.6  E-5,  and  cx  =  1 

using  sixth  order  polynomials  to  fit  the  parameter  histories. 
The  converged  parameter  histories  are  shown  in  Figures  5,  6, 
and  7  along  with  the  initial  parameter  histories.  These 
figures  display  the  effectiveness  of  the  method  used  to 
calculate  the  initial  parameters .  The  curves  show  corre¬ 
sponding  peaks  and  valleys  and  are  of  the  same  order,  which 
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Table  I 


State  Histories  (Run  l) 


Time 
( sec ) 

x,  of  G 
x,  of  Y 
x,  -DIFF 
(ft) 

xa  of  G 
xa  of  Y 

Xj -DIFF 
(ft) 

x3  of  G 
x}  of  Y 
X3-DIFF 
(ft) 

xv  of  G 
x,  of  Y 
x» -DIFi 
(ft/sec ) 

xr  of  G 
Xy  Of  Y 
xr-DIFF 
(rad) 

x4  of  G 
x4  of  Y 
x* -DIFF 
(rad) 

6.0 

2515^-3 
;  25114.3 

-40. 

12182.7 

12285.1 

102.4 

8049.5 

8085.3 

36.8 

582.2 

575.7 

-6.5 

4.26021 

4.2481 

-.01211 

-.352506 

-.34208 

.010526 

18.0 

19323.5 

19219.4 

-104.1 

9539.4 

9652.1 

112.7 

5724.3 

5688.8 

-35*5 

552.9 

562. 

9.1 

4.29792 
4.2603 
-.03762  j 

-.345674 
-.35954 
-.01 3666 

30.0 

1  3237 . 8 
1341  5.6 
-322.2 

6967.9  1 
6982.5 
14.6  ■ 

3504.7 

3533.5 

28.8 

541.5 

531.6 

-9.7 

4.25736 
4.2324 
-.0249 6 

-.343423 

-.32987 

.013553 

42.0 

8231  .2 
7Q7P .  5 
-252.7 

4137.8 

4112.8 

-25. 

1482. 

1336.6 

-145.4 

542.6 

543.6 

1 . 

4.22642 

4.281 

-.02542 

-.269352 
-.33685 
-.06749 : 

54.0 

2710.5 

2400.7 

-309.8 

1279.8 

1 1 p  2 . 3 

-97.5 

282.5 

142.2 

-140.3 

502.7 
500 . 6 

-2.1 

4.25053 

4.2726 

.02207 

-.102-91 

-.04538 

.057111 

66 . 0 

-2228.9 
-2527.5 
-2QP .  6 

-1074. 
-ll^C. 2 
-105.2 

74.5 

45.4 

-29.1 

405.4 

402.7 

-2.7 

4.2807 
4.26  ;« 
-.0169 

.003491 

-.01097 

-.013561 

78.0 

-6l  42 , 
-6448.9 
-305.9 

-2892.5 

-3066.7 

-174.2 

55.3 

14.9 

-40.4 

321  .8 
320.8 

-1 . 

4.26536 

4.2726 

.00724 

1  - 

- . 0331 36 
-.00524 
.002898 

DIFF  =  Difference  between  of  Y  and  xu  of  G 
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Fig.  5.  Initial  and  Converged  CL  (Hun  l) 


itial  and  Converged  (Run  l) 


**# 


Initial  and  Converged  /-*  (Run  l) 


1 


means  the  initial  parameter  histories  did  not  require  exten¬ 
sive  changes,  implying  fast  convergence.  Although  the 
magnitudes  differ  somewhat  for  CL  and  CD,  this  can  be 
attributed  to  the  inaccuracy  of  the  reference  area  s(r^kou~‘)  and 
the  average  value  used  for  w  (198,200  lb),  which  were  used 
in  calculating  the  initial  parameters.  The  accuracy  of  the 
algorithm  for  this  run  is  on  the  order  of  100  ft  for  the 
xi,  X2 /  and  X3  states,  5  ft/sec  for  X4,  and  0.02  rad  for 
states  X5  and  x^  (Table  I)  . 

Run  2 

The  second  run  spanned  the  first  thirty-nine  seconds 
of  Run  1  and  used  500  CP  seconds  to  converge  to  A  =  5.2  E-9, 

B  =  7.9  E-9,  and  ex  =  1 .  The  order  of  the  polynomial 

used  to  fit  the  parameter  histories  was  varied  from  second 
to  sixth  order.  Figure  8  shows  the  converged  lift  coeffi¬ 
cient  histories  for  the  two  extreme  cases  of  second  order 
and  sixth  order.  The  second  order  curve  slices  through 
the  hills  and  valleys  of  the  sixth  order  curve,  yet,  as 
Tables  II  and  III  show,  the  state  history  fits  for  both 
runs  are  nearly  the  same.  This  illustrates  the  importance 
of  choosing  a  polynomial  of  large  enough  order  such  that 
all  of  the  gross  characteristics  of  the  true  parameter 
history  are  mapped. 

This  interval  was  also  used  as  a  comparison  to  Run  1 
to  study  how  choice  of  interval  length  and  order  of  poly- 
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Second  Order  Polynomials  (Run  2) 


Table  II 


Second  Order  Polynomial  State  Histories  (Run  2) 


1  5.0 


21  .0 


27.0 


33.0 


39.0 


x,  of  G 
x,  of  Y 
x, -DIFF 
(ft) 

x*  of  G 
Xj  of  Y 
x^-DIFF 
(ft) 

x,  of  G 
x*  of  Y 
x, -DIFF 
(ft) 

x*  of  G 
xi  of  Y 
x* -DIF  F 
(ft/sec ) 

26636. 5 

12931 .6 

8648.4 

580.1 

26603.4 

12965.6 

6667.5 

582.5 

-33.1 

34. 

19.1 

2.4 

2  37  72.2 

11551. 

7400.4 

570.4 

23532.4 

11621.9 

7504.6 

565.3 

-rQ.8 

70.9 

104.2 

-5.1 

20834.5 

10240.7 

6204.2 

562.1 

20503.1 

10312.7 

6291. 

560.4 

-141.4 

72. 

86.8 

-1.7 

17°5R .4 

8978.8 

5100.8 

552.5 

17747.8 

8986.9 

5108.9 

558.9 

-210.6 

48.1 

8.1 

6.4 

15137.1 

7606.9 

4075.3 

542.4 

14830.2 

7656. 

4045.3 

540.5 

-306.9 

49.1 

-30. 

-1.9 

12398.9 

6227. p 

3061  .4 

53^*8 

12038.7 

6282.2 

3009.4 

529.6 

-359.2 

54.4 

-52. 

-5.2 

9754.7 

4799.2 

1957.1 

534. 

9328.2 

4862.2 

1899.9 

538.8 

-426.5 

62.3 

-57.2 

_ 

4.8 

xf  of  G 
x,  of  Y 
Xc-DIFF 
(rad) 


4.25241 

4.2394 

-.01301 

4.2516 

4.2569 

-.024? 

4.2-054 


4.28144 


4.23146 


012618 


4.2595?  I -.313999 


012066 


022719 


Difference  between  xt  of  Y  and  of  G 


Table  III 


Sixth  Order  Polynomial  State  Histories  (Run  2) 


Time 

(sec) 

x  of  G 
x  of  Y 
x  -DIF? 
(ft) 

i - 

x  of  G 
x  of  Y 
x  -DIFF 
(ft) 

x  of  G 
x  of  Y 
x  -DIFF 
(ft) 

x  of  G 
x  of  Y 
x  -DIFF 
(ft/sec) 

x  of  G 
x  of  Y 
x  -DIFF 
(rad) 

x  of  G 
x  of  Y 
x  -DIFF 
(rad) 

3.0 

266  76. 5 
26607. 4 
-33.1 

12931 .6 

12965.6 
34. 

8648.4 

8667.5 
19.1 

580.1 

582.5 

2.4 

4.25241 

4.2394 

-.01301 

- . 366667 

-.34208 

.024587 

9.0 

23732.2 

28603.4 

-OQ,P 

11  851  . 
11621.9 
70.9 

740C.4 

7504.5 

104.2 

570.4 

565.3 

-5.1 

4.2815 

4.2569 

-.0247 

-.358668 

-.35605 

.012618 

15.0 

20834. 5 
20693.1 
-1 4l  .4 

10240.7 

10312.7 
72. 

6204.2 
6291  . 
86.8 

562.1 

560.4 

-1.7 

4.29054 

4.2638 

-.02674 

-.346619 

-.37001 

-.021391 

21.0 

17g58.4 

17747.8 

-210.6 

8938.8 

COP5.9 

48.1 

5100.8 

5108.9 

8.1 

552.5 

558.9 

6.4 

4.28144 

4.2534 

-.02804 

-.324627 
-.34208 
-. 01^453 

27.0 

15132.1 

14830.2 
-306.9 

7605.9 

7656. 

49.1 

4075-3 

4045.3 

-30. 

542.4 

540.5 

-1.9 

4.35957 

4.2603 

.0C073 

-.314 

-.32114 

-.00714 

33.0 

12398.9 

12039.7 

-359.2 

6227.8 

6282.2 

54.4 

3061  .4 
3009.4 

-52. 

534.8 

529.6 

-5.2 

|  4.23146 
:  4.2185 
-.01296 

-.330014 

-.34208 

-.012066 

39.0 

9754.7 

0328.2 

-426.5 

4799.9 

4862.2 

62.3 

1957.1 
1899. 9 
-57.2 

534. 

538.8 

4.8 

4.2035 

4.2045 

.001 

-.380419 

-.3577 

.022719 

DIFF  =  Difference  between  x;.  of  Y  and  of  G 
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nomial  affects  the  accuracy  and  convergence  of  the  algorithm. 
A  comparison  of  Tables  I  and  II,  both  calculated  using 
sixth  order  polynomials,  shows  that  Run  2  was  only  slightly 
more  accurate  than  Run  1  at  fitting  the  experimental  state 
histories  over  the  same  time  intervals.  Yet,  comparing 
Figures  9,  10,  and  11  to  the  Run  1  figures  shows  that  the 
true  parameter  histories  are  more  fully  described  by  Run  2. 
Therefore,  the  appropriate  length  of  the  interval  and  the 
polynomial  order  are  very  dependent  upon  each  other  for 
properly  describing  the  true  parameter  histories.  Again, 
the  initial  parameter  histories  give  a  good  means  for 
choosing  an  appropriate  polynomial  order. 

This  interval  was  also  used  to  test  the  effect  of 
different  W  diagonal  elements.  Tables  IV-VI  represent 
three  state  histories,  each  with  a  different  W.  Table  IV 
uses  1.E4  for  the  first  three  diagonal  elements  of  W,  (call 
them  wl,  w2,  w3,  respectively),  1.  for  w4 ,  l.E-4  for  w5 
and  w6 .  This  represents  very  accurate  radar  measurements 
for  the  x^ ,  x2,  and  x^  states  relative  to  the  measurements 
for  the  x5  and  xg  states.  Table  V,  W2 ,  is  with  all  six 
diagonal  elements  as  1,  such  that  all  six  radar  measurements 
are  assummed  of  equal  accuracy.  Table  III  is  the  opposite 
extreme  to  Table  IV  with,  wl ,  w2,  w3  as  l.E-8,  w4  as  l.E-4, 
w5  and  w6  as  1.  Comparing  the  convergence  characteristics 
of  the  algorithm  using  these  different  weighting  matrices 


Fig..  10.  Initial  and  Converged  Cp  (Run  2) 


TIME  HISTOR 


3: 


Fig.  It .  Initial  and  Converged^  (Hun  2) 


Table  IV 

State  Histories  for  Wl  (Run  2) 


Time 

(sec) 

x,  of  0 
x,  of  Y 
x.  -DIFF 
(ft) 

x*  of  G 
xa  of  Y 
Xa-DIFF 
(ft) 

x5  of  G 

X,  of  Y 
x.-DIFF 
(ft) 

x.  of  G 
xH  of  Y 
xij-DIFF 
(ft/sec) 

x,  of  G 
of  Y 
X5-DIFF 
(rad) 

xk  of  G 
x*  of  Y 
x* -DIFF 
(rad) 

3.0 

26607.2 

26603.4 

-3.9 

I203fc.7 

12065.6 

25.9 

8678.1 

8667.5 

-10.6 

584.9 

582.5 

-2.4 

4.27907 

4.2394 

-.03967 

-.330665 

-.34208 

-.011415 

9.0 

23606.9 

2?632.4 

25.5 

11641.5 
11621 .9 
-10.6 

7553.5 

7504.5 
-48.9 

569.6 

565.3 

-4.3 

4.30936 

4.2569 

-.05246 

-.344296 

-.35605 

-.01175- 

15.0 

20683.5 

20693.1 

9.5 

10348. 

10312.7 

-35.3 

6331.5 

6201. 

-40.5 

573.8 

560.4 

-13-4 

4.28289 

4.2638 

-.01909 

-.376847 

-.37001 

.008837 

21  .° 

17^53.2 

17747.8 

-5-4 

8094.8 

8985.9 

-7.9 

5089.7 

51 08. 9 
19.2 

575.7 

558.9 

-16.8 

4.2808 

4.2534 

-.0274 

-.343946 

-.34208 

.001866 

27.0 

14022. 3 
14830.2 
7.9 

7636.7 

7655. 

19.3 

4032.2 

4045.3 
13.1 

554.2 

540.5 

-13.7 

4.27064 

4.2603 

-.01034 

-.29=°67 

-.32114 

-.022273 

33.0 

12044.3 

12039.7 

-4 . 6 

6283. 

6282.2 

-.8 

3017.5 

3009.4 

-8.1 

534.3 

529.6 

-4.7 

4.24955 

4.2185 

-.03105 

-.350458 

-.34208 

.006378 

39.0 

9368.6 

9329.2 

-40.4 

4848 . 3 
4862.2 
13.9 

1915. 1 
1899. 
-15-2 

548.4 

538.8 

-9.6 

4.1116 

4.2045 

.0929 

-.283106 

-.3577 

-.074594 

DIFF  =  Difference  between  xj_  of  Y  and  of  G 
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Table  V 

State  Histories  for  W2  (Run  2) 


Time 

(sec) 

x,  of  G 
x,  of  Y 

X|  -DIFF 
(ft) 

Xi  of  G 
x*  of  Y 
x^-DIFF 
(ft) 

x*  of  G 

Xi  of  Y 

Xi -DIFF 
(ft) 

x*  of  G 
xi  of  Y 
x#-DIFF 
(ft/sec) 

x{  of  G 
x*  of  Y 
x- -DIFF 
(rad) 

xfc  of  G 
xu  of  Y 

Xt -DIFF 
(rad) 

3-0 

265P0.7 

26603.4 

22.7 

12025*6 

12965.6 

40. 

8658.1 

8667.5 

9.4 

5n7.9 

582.5 

-15.4- 

4.27551 

4.2394 

-.03611 

-.338229 

-.34208 

-.003851 

9.0 

23604.9 

23632.4 

27.5 

11615.6 

11621.9 

6.3 

7509.5 

7504.6 

-4.9 

558.1 

565-3 

7.2 

4.30572 

4.2569 

-.04882 

-.350568 

-.35605 

-.005462 

15.0 

20716.3 

20603.1 

-23.2 

10343.6 
10  712.7 
-30.9 

6309.9 
6291  • 
-18.9  j 

573.8 

560.4 

-13.4 

4.28822 

4.2638 

-.02442 

-.357956 

-.37001 

-.702052 

21.0 

17^41.3 
17747. 8 

6.5 

P985.2 

8085.9 

1.7 

5102.9 

5108.9 

6. 

580. 

558.9 

-21.1 

4.28332 

4.3534 

-.02Q92 

-.333341 

-.34208 

-.008739 

27.0 

14819.2  | 

14830.2 
11 . 

764-7.1  1 
7655. 
11.9 

4045.5 

4045.3 

-.2 

1  546.6 

5^0.5 

i  -6.1 

4.27632 

4.2603 

-.01602 

-.3112-2 

-.32114 

-.009858 

33.0 

12056.2 

12039-7 

-16.5 

62-8.9 

6282.2 

-16.7 

3010.7 

3009.4 

-1.3 

544.1 

529.6 

-14.5 

4.23853 

4.2185 

-.02003 

-.342631 

-.34208 

-.000551 

39.0 

9325.3 
032°. 2 
2.9 

4848.9 
4862 .2 
13.3 

1905.2 

18°9.9 

-6.3 

522.9 

538.8 

15.9 

4.22591 

4.2"45 

-.02141 

-.322315 

-.3577 

-.035365 

DIFF  =  Difference  between  of  Y  and  of  G 
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Table  VI 


State  Histories  (Run  3) 


Time 

(sec) 

x,  of  G 
x,  of  Y 
x,  -DIFF 
(ft) 

x,  of  G 
Xi  of  Y 
x  j -DIFF 
(ft) 

1 - 

x  j  of  G 

Xj  of  Y 
x,-DIFF 
(ft) 

x,  of  G 
x,  of  Y 
x^ -DIFF 
(ft/sec) 

xr  of  G 
Xj-  Of  Y 

xs -DIFF 
(rad) 

xL  of  G 
xt  of  Y 
xL -DIFF 
(rad) 

3.0 

-53292.5 

-6l67.2 

17065.8 

605.7 

I 

-1 .19648 

-.255783 

-53262.6 

-6175.3 

17064.7 

613.1 

-1.21249 

-.250091 

29.8 

-8.1 

-1.1 

7.4 

-.016 

.006 

12.0 

-49293.1 

-2820.2 

15593.5 

605.2 

-1 .46558 

-.297804 

-49092.5 

-2987.3 

15695.2 

599.2 

-1 .47314 

-.278401 

200.6 

-167. 

101.7 

-5.9 

-.008 

.019 

21  .0 

-44^26.7 

-856.4 

13870.4 

606.2 

-1.76756 

-.341253 

-44*08.9 

-872.7 

13913.3 

609.9 

-1.67224 

-.381764 

«?.e 

-16.3 

42.8 

3.7 

.095 

-.041 

30.0 

-39500.2 

-32.O 

11957.9 

603.5 

-1 .91483 

-.388649 

-39482.7 

102.3 

H881.5 

602.9 

-1 .09473 

-.381707 

17.4 

135-1 

-76. 3 

-.6 

.01 

-.004 

39.0 

-3451  3.3 

209.2 

9915.7 

594.5 

-1.98017 

-.365888 

-34461.1 

2  30.7 

9  Q99. 

592.1 

-1.99489 

-.359206 

52.1 

-68.5 

82.3 

-2.4 

-.015 

.017 

48.0 

-20584.3 

375.4 

8053.8 

576.2 

-2.00653 

-.332453 

-29512.6 

256.9 

8217.4 

578. 

-1.97712 

-.351077 

71  .7 

-118.5 

158.5 

1.8 

.029 

-.019 

57.0 

-24736.1 

510.8 

6301.3 

575-5 

-1.98293 

-.410185 

-24670. 

3P  5.1 

6461 . 3 

573-2 

-1.9584 

-.376647 

66 . 1 

-125.7 

160.4 

-2.3 

-.002 

.039 

DIFF  =  Difference  between  of  Y  and  x\_  of  G 
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provides  a  method  for  choosing  the  best  weighting  matrix. 

As  would  be  expected,  the  heavily  weighted  states 
converged  better  for  all  three  cases.  Table  V  shows  that 


the  algorithm  ignores  the  small  x^  and  Xg  histories 
because  they  contribute  little  to  the  performance  index 
compared  to  the  x^,  X2 ,  and  X3  states.  Table  III  shows 
the  order  of  weighting  according  to  the  Appendix,  where 
the  range  (x^,  X2  and  X3)  has  uncertainties  on  the  order 
of  50  ft  and  the  azimuth  and  elevation  measurement  (X5  and 
xg)  uncertainties  are  five  to  seven  magnitudes  smaller  at 
2.E-3  rad  to  2.E-5  rad. 

Table  IV  has  the  best  fit  for  the  states  x^,  x2 ,  and 
X3  but  the  algorithm  never  converged.  This  can  be  explained 
by  the  extremely  small  weights  for  x5  and  xfe ,  thereby 
causing  the  algorithm  to  ignore  these  states  when  computing 
the  performance  index.  For  this  reason,  the  parameter 
histories,  which  are  very  dependent  upon  the  X5  and  x6  states 
(see  the  state  equations) ,  did  not  converge  to  the  proper 
values.  The  same  reasoning  applied  to  the  Table  V  weights, 
for  which  the  run  did  not  converge.  Therefore,  although 
the  weighting  matrix  for  Table  III  had  a  less  accurate  fit 
to  the  experimental  trajectory  states  x^,  X2/  and  x^,  it  is 
the  matrix  that  most  closely  represents  the  true  error  model. 
The  errors  for  the  states  in  Table  III  are  representative  of 
the  magnitude  of  the  errors  that  can  be  achieved  by  this 


algorithm  while  providing  converged  (good)  parameter 
histories . 

Run  3 

The  third  run  used  data  from  the  BET  as  the  experi¬ 
mental  values.  As  was  previously  mentioned,  the  BET 
provides  state  history  data  which  was  used  in  this  run. 

Since  the  BET  also  provides  the  CL  and  CD  histories  for 
this  trajectory  data,  the  converged  algorithm  should  yield 
similar  parameter  histories.  This  allows  the  algorithm  to 
"prove"  itself  against  the  BET's  parameter  estimation 
procedure . 

This  run  used  500  CP  seconds  to  converge  to  A  =  l.E-7, 

B  =  6.E-5,  and  o<  =1.  The  state  time  histories  are 
listed  in  Table  VI.  The  initial  and  converged  parameter 
histories  are  portrayed  in  Figures  12,  13,  and  14.  Notice 
that  the  bank  angle  fit  was  very  good,  as  the  accuracy  of 
the  Xj.  state  history  fit  implies,  although  the  initial  and 
converged  parameter  histories  are  significantly  different 
in  magnitude.  This  shows  the  ability  of  the  algorithm  to 
converge  given  a  poor  initial  yS  . 

Because  the  lift  to  drag  ratio  (L/D)  can  be  calculated 
without  the  need  for  a  reference  area,  s,  and  a  weight,  w, 
the  L/D  for  the  converged  f3  was  compared  to  the  L/D  values 
given  by  NASA's  BET,  so  that  the  comparison  is  only  dependent 
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Initial  and  Converged  CD  (Run  3) 


V  pa 


upon  the  two  methods  of  parameter  estimation.  Figures  15 
and  16  present  in  graphic  form  the  comparisons  of  the 
initial  to  the  converged  and  the  converged  to  the  BET  L/D 
ratios,  respectively.  Figure  16  shows  good  correspondence 
between  the  two  methods  of  identifying  the  L/D  ratio.  The 
BET  L/D  history  is  not  restricted  to  be  a  smooth  polynomial. 
The  beginning  and  end  of  the  converged  polynomial  L/D 
history  does  not  correspond  as  well  to  the  BET  values. 

This  can  be  attributed  to  the  assumption  that  xQ  was  perfectly 
known  and  to  the  use  of  polynomials  of  lower  order  than 
may  be  necessary  for  the  interval  length  (Figure  12,  for 
instance) . 
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Fig.  15.  Initial  and  Converged  L/D  (Run  3) 


V  Conclusion 


Summary 

This  report  presented  a  nonlinear  parameter  identi¬ 
fication  algorithm  that  was  used  to  find  C^,  CD,  and 
time  histories,  derived  from  radar  trajectory  data,  for 
the  Space  Shuttle  reentry.  The  algorithm,  a  Square-Root 
Variable-Metric  optimization  technique,  (Ref  4,9)  provided 
a  fast  and  inexpensive  method  for  obtaining  good  results 
for  the  parameter  histories,  by  minimizing  the  weighted 
least  squares  performance  index,  J,  of  the  experimental 
versus  calculated  state  histories. 

There  were  several  programming  considerations 
discussed  in  Chapter  III.  The  procedure  used  to  execute 
the  one  dimensional  search  was  detailed,  showing  a  need 
to  manipulate  the  tolerances  for  convergence  and  to  vary 
the  step  size  when  needed.  Finding  the  best  numerical 
differentiation  constant,  €.  ,  was  done  by  varying  € 

until  the  most  accurate  estimate  of  the  derivative  was 
obtained.  This  reduced  errors  due  to  computer  round-off 
and  central  difference  formula  truncation.  Nondimensional 
polynomial  coefficients  were  also  used  to  reduce  these 
errors.  Finally,  by  inverting  the  magnitude  of  the 
uncertainties  given  in  the  Appendix,  the  diagonal  elements 
of  the  W  matrix  were  chosen. 

The  results  were  presented  and  analyzed  in  Chapter  IV. 


Run  1  and  Run  2  tested  the  mechanics  of  the  algorithm 
dependent  upon  the  choice  of  interval  length,  polynomial 
order,  and  weighting  matrix.  It  was  shown  that  the  initial 
parameter  histories  were  close  enough  to  the  converged 
values  to  give  fairly  rapid  convergence.  Run  3  used  the 
BET  parameters  as  a  comparison  to  the  converged  solution 
using  the  BET  state  histories  as  the  radar  trajectory  data. 

This  showed  good  correspondence  between  the  two  estimation 
methods  for  a  L/D  ratio  time  history,  hence,  giving  confi¬ 
dence  in  the  algorithm's  effectiveness  (Fig.  16).  It  was 
noted  that,  because  of  the  characteristics  of  polynomial 
approximation  and  the  assumption  that  xQ  was  fully  known, 
the  converged  parameter  histories  may  be  inaccurate  at 
the  ends  of  the  time  intervals.  All  three  runs  converged 
well  to  the  given  experimental  state  histories  (Tables  I, 

II,  and  III) ,  but  not  all  of  the  parameter  history  charact¬ 
eristics  were  portrayed  if  the  order  of  the  polynomial  fit 
did  not  match  correctly  with  the  interval  length. 

Computer  time  was  short,  about  500  CP  seconds  (hence, 
inexpensive)  and  was  made  shorter  by  choosing  smaller  order 
polynomials . 

Recommendations 

This  algorithm  provides  a  relatively  simple  method  for 
identifying  parameter  histories.  The  results  for  the  nonlinear 
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state  model,  deterministic,  Square-Root  Variable-Metric 
optimization  routine  corresponded  well  with  the  BET  results 
without  the  problems  of  matrix  inversion.  The  computer 
load  is  very  manageable  with  500  CP  seconds  as  an  average 
time  for  fitting  about  a  minute  long  interval  with  a  sixth 
order  polynomial  using  the  CDC  Cyber  74. 

Because  the  radar  data  is  readily  available  for  most 
of  the  Space  Shuttle  reentry  trajectory,  this  algorithm 
could  be  useful  for  quickly  yielding  good  parameter  history 
results.  These  results  would  allow  evaluation  of  maneuvers 
or  unplanned  motion  without  the  long  wait  for  NASA's  BET. 
This  could  perhaps  be  useful  for  speeding  up  the  turn 
around  time  between  launches  should  any  problems  which 
might  arise  during  reentry  need  evaluation  before  the 
next  launch. 

Redundant  results  are  important  when  testing  the 
performance  of  a  system.  Using  this  parameter  identifica¬ 
tion  process  would  give  a  second,  redundant,  means  for 
evaluating  the  accuracy  of  the  BET  data.  Although  NASA's 
BET  gives  a  full  six  DOF  evaluation  of  the  Space  Shuttle’s 
reentry,  such  a  comparison  would  inspire  more  confidence  in 
the  BET  results. 

The  method  used  here  to  identify  parameters  can  also 
be  used  where  any  parameters  need  to  be  found,  as  long  as 
an  appropriate  state  model  is  given.  The  polynomial  fits. 
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one-dimensional  search  development,  and  weighting  matrix 
choice  are  all  appropriate  considerations  whether  using 
this  or  a  similar  optimization  technique.  Therefore,  the 
results  presented  give  insight  to  what  can  be  achieved 
in  the  field  of  nonstochastic  parameter  identification. 

Further  work  with  this  algorithm  could  explore  the 
use  of  a  full  measurement  covariance  matrix  inversion 
as  the  weighting  matrix.  This  would  more  accurately 
describe  the  real  world  radar  errors  and  thereby  improve 
the  accuracy  of  the  algorithm.  Also,  the  state  model 
could  be  expanded  to  more  accurately  describe  the  Space 
Shuttle  three  DOF  motion.  Spherical  Earth,  wind,  and 
noise  modelling  would  all  increase  the  accuracy  of  the 
model,  but  at  some  degree  of  expense  to  the  speed  of 
convergence,  as  the  computational  load  increases 
due  to  a  more  complex  dynamical  system  model.  The  trade-off 
of  accuracy  versus  speed  of  convergence  can  be  studied 
using  the  results  of  these  further  studies  and  the  results 
presented  in  this  report. 
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Appendix 


The  following  four  tables  are  from  Reference  8. 

The  first  two  tables,  Cr%  and  C“2. »  show  the  weights  and 
apriori  bias  information  for  the  Multisensor  Report  data. 
Tables  D-l  and  DJshow  the  radar  measurement  mean  and 
uncertainties.  The  diagonal  elements  for  the  weighting 
matrix  were  arrived  at  by  combining  the  bias  and  uncertainty 
magnitudes  and  inverting  them.  The  first  three  diagonal 
elements,  representing  states  x^ ,  X2,  and  X3,  are  weighted 
from  the  range  magnitude  and  the  last  two  diagonal  elements, 
representing  x5  and  xg,  are  weighted  from  the  azimuth  and 
elevation  magnitudes.  The  fourth  diagonal  element,  repre¬ 
senting  the  weighting  of  the  xA  (velocity)  state,  was 
given  a  value  of  magnitude  between  the  range  and  azimuth 
magnitudes  since  Doppler  residual  statistics  were  not 
provided. 
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SENSOR  WEIGHTS  AND  APRIORI  BIAS  INFORMATION  FOR  RUN  NO. 


WEIGHTS 

BIASES 

BIAS  UNCERTAINTIES 

SENSOR 

RANGE 

(FEET) 

AZIMUTH 

(DEG) 

ELEV. 

(DEG) 

RANGE 

(FEET) 

AZIMUTH 

(DEG) 

ELEV. 

(DEG) 

RANGE 

(FEET) 

AZIMUTH 

(DEG) 

ELEV. 

.  (DEG) 

FPS16/38 

30.00 

.008 

.008 

19.52 

.00474 

.00837 

63.09 

.03261 

.03266 

NASA 

FPS16 

30.00 

.008 

.008 

10.41 

.00065 

.01140 

58.62 

.03307 

.02606 

TPQ18 

SA  48 

30.00 

.008 

.008 

-38.25 

.01682 

.01991 

65.66 

.02749 

.02993 

Cl 

0.00 

.003 

.003 

0.00 

-.00143 

-.00163 

0.00 

.00524 

.00344 

C4 

0.00 

.003 

.003 

0.00 

.00318 

-.00159 

0.00 

.00407 

.00383 

C5  • 

0.00 

.003 

.003 

0.00 

-.00139 

-.00197 

0.00 

.00496 

.00287 

C6 

0.00 

.003 

.003 

0.00 

.00574 

-.00191 

0.00 

.01010 

.00435 

C7 

0.00 

.003 

.003 

0.00 

-.00379 

.00296 

0.00 

.00627 

.00509 

TABLE  C-l 


BIAS  ESTIMATION  RESULTS  AND  UNCERTAINTIES  FOR  RUN  NO.  1 

SENSOR 

BIASES 

BIAS  UNCERTAINTIES 

JtANGE 

(FEET) 

AZIMUTH 

(DEG) 

ELEVATION 

(DEG) 

RANGE 

(FEET) 

AZIMUTH 

(DEG) 

ELEVATION 

(DEG) 

FPS16/38 

27.87364 

.00109 

.02611 

1.38904 

. 00035 

. 00039 

NASA  FPS16 

17.88180 

-.00161 

.02940 

1.29678 

.00033 

.00039 

TPQ18  SA  48 

-27.56410 

.01939 

.04245 

1.76858 

-00038 

.00044 

Cl 

0.00000 

-.00712 

. 00680 

.00000 

.00026 

.00033 

C4 

0.00000 

.00281 

.00476 

.00000 

.00019 

.00026 

C5 

0.00000 

-.00280 

.00504 

.00000 

.00021 

.00027 

C6 

0.00000 

.00632 

.00498 

.00000 

.00020 

.00025 

C7 

0.00000 

-.01127 

.01350 

.00000 

.00029 

.00036 

TABLE  C-2 


RESIDUAL  STATISTICS  FOR  RUN  NO.  1 


SENSOR 

MEAN 

UNCERTAINTIES 

RANGE 

(FEET) 

AZIMUTH 

(DEG) 

ELEVATION 

(DEG) 

RANGE 

(FEET) 

AZIMUTH 

(DEG) 

ELEVATION 

(DEG) 

FPS16/38 

.00273 

-.00000 

.00000 

61.20924 

.03209 

.02937 

NASA  FPS16 

.00252 

-.00000 

.00000 

59.67920 

.03065 

.02285 

TPQ18  SA  48 

.00404 

.00000 

.00000 

64.70393 

.02703 

.02392 

Cl 

- 

-.00000 

.00002 

- 

.00521 

.00315 

C4 

- 

-.00000 

.00001 

- 

.00438 

.00392 

C5 

-.00000 

.00002 

- 

.00476 

.00285 

C6 

- 

.00000 

.00001 

- 

.01114 

.00439 

C7 

- 

-.00000 

.00001 

- 

.00639 

.00427 

TABLE  D-l 


POSITIONAL  UNCERTAINTY  STATISTICS  FOR  RUN  N 

0.  1 

COORDINATE 

MEAN  (FEET) 

UNCERTAINTIES  (FEET) 

R  SPHERICAL  ERROR 

43.76213 

60.96327 

X  TOTAL 

'  14.33824 

15.27715 

Y  TOTAL 

25.19341 

38.38616 

Z  TOTAL 

31.88084 

45.47689 

X  RANOOM 

14.31515 

15.25862 

Y  RANDOM 

25.15378 

38.33737 

Z  RANOOM 

31.80159 

45.39780 

X  SYSTEMATIC 

.80540 

.76054 

Y  SYSTEMATIC 

1.39732 

1.94573 

Z  SYSTEMATIC 

2.17653 

2.73815 

TABLE  D-2 
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